%
% Program to solve transmission line equations
%
% by Joe LoVetri, 2010
%

% set up physical parameters
%
close all
clear all

c_0 = 3e8;          % speed of propagation
C = 10e-12;         % PUL capacitance
L = 1/c_0^2/C;      % calculate PIL inductance based on speed and C
Len = 100;          % length of line

Rt_near = 0;
Vt_near = 0;
Pulse_begin = 10e-9;
Pulse_end = 200e-9;

Rt_far = sqrt(L/C);
Vt_far = 0;

% set up numerical parameters

CFL = 1.0;                  % CFL stability factor
num = 100 ;                  % number of space points
delta_x = Len/num;          % space step size
delta_t = CFL*delta_x/c_0;  % time step based on CFL stability condition.
xv = [0:num]*delta_x;             % vector of space locations for voltage
xi = [1:num]*delta_x-delta_x/2;   % vector of space locations for current

% initialize voltage array and current array

current = zeros(1, num);
voltage = zeros(1, num+1);

% voltage(45:55) = 1;                    % square pulse as initial condition
%voltage(45:55) = [(0:1:5) (4:-1:0)];   % square pulse as initial condition
% Gaussian distributed initial voltage
x=1:num+1;
voltage = exp(-((x-50).^2)/50);

% plot the initial voltage and current
subplot(2,1,1)
plot(xi, current)
xlim([min(xi), max(xi)]);
xlabel('Length along Tx-Line [m]')
ylabel('Current')
title(['Time = ' num2str(0)  ' [nsec]']);

subplot(2,1,2)
plot(xv, voltage)
xlabel('Length along Tx-Line [m]')
ylabel('Voltage')
title(['Time = ' num2str(0) ' [nsec]']);

reply = input('More? y/n [y]: ', 's');
    if reply == 'n' 
        return 
    end;    

n = 1;
time_v = n* delta_t;
time_i = time_v - delta_t/2;
factor = delta_t/delta_x;

% first update is different for the current (need to divide factor by 2)

% Boundary conditions
voltage(1) = Vt_near - current(1)*Rt_near;
voltage(num+1) = Vt_far + current(num)*Rt_far;

% Updat interior points
current(1:num) = current(1:num) - ...
    (factor/L)*(voltage(2:num+1) - voltage(1:num));
voltage(2:num) = voltage(2:num) - ...
    (factor/C)*(current(2:num) - current(1:num-1));


% plot the first update of voltage and current
subplot(2,1,1)
plot(xi, current);
xlim([min(xv), max(xv)]);
xlabel('Length along Tx-Line [m]');
ylabel('Current');
title(['Time = ' num2str(round(time_i*1e10)/10) '[nsec]']);

subplot(2,1,2);
plot(xv, voltage);
xlabel('Length along Tx-Line [m]')
ylabel('Voltage')
title(['Time = ' num2str(round(time_v*1e10)/10) '[nsec]']);

reply = input('More? y/n [y]: ', 's');
    if reply == 'n' 
        return 
    end;    
    
stop = 0;
while stop == 0 
    n = n+1;
    time_v = n* delta_t;
    time_i = time_v - delta_t/2;

    % Boundary Conditions
    % calculate Thevenin voltage
    if time_v > Pulse_begin && time_v < Pulse_end
        Vt_near = 0;
    else
        Vt_near = 0;
    end
    voltage(1) = Vt_near - current(1)*Rt_near;
    voltage(num+1) = Vt_far + current(num)*Rt_far;
    
    % Update interior points
    current(1:num) = current(1:num) - ...
        (factor/L)*(voltage(2:num+1) - voltage(1:num));
    voltage(2:num) = voltage(2:num) - ...
        (factor/C)*(current(2:num) - current(1:num-1));
    
    % Plot
    subplot(2,1,1)
    plot(xi, current);
    xlim([min(xv), max(xv)]);
    xlabel('Length along Tx-Line [m]');
    ylabel('Current');
    title(['Time = ' num2str(round(time_i*1e10)/10) ' [nsec]']);
    
    subplot(2,1,2);
    plot(xv, voltage);
    xlabel('Length along Tx-Line [m]')
    ylabel('Voltage')
    title(['Time = ' num2str(round(time_v*1e10)/10) ' [nsec]']);
    
    reply = input('More? y/n [y]: ', 's');
    if reply == 'n' 
        stop = 1 
    end;    
end 